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We develop a general theory for thermal transport in anharmonic systems under the weak system- 
bath coupling approximation similar to the quantum master equation formalism. A current operator 
is derived, which is valid not only in the steady state, but in the transient regime as well. Here 
we focus on the effects of anharmonicity on the steady state thermal conductance of a mono- and 
di-atomic molecular system. We also study molecules being confined in a double well potential. 
We find that when the molecules have a non-linear on-site potential the low-temperature thermal 
conductance is dramatically affected by the strength of non-linearity, whereas for the diatomic 
molecule connected by an anharmonic spring the strength of anharmonicity plays almost no role in 
the low-temperature regime. In case of the molecules confined in a double well potential we find 
that the height of the barrier greatly affects the thermal conductance; once the molecules can feel 
the effect of the barrier we observe negative differential thermal conductance at both high and low 
temperatures. 



I. INTRODUCTION 

The theory of thermal transport dates back to the 
works of Debye and Peierlsi^ who studied the heat trans- 
fer within solids. In recent years the study of heat trans- 
fer in nano systems is very active due to the need in 
device applications. One approach is to look at the 
classical properties of heat transport using molecular 
dynamics^. This technique gives good insight in the high- 
temperature regime but cannot be applied to low temper- 
atures. The low-temperature regime can be probed us- 
ing noncquilibrium Green's function (NEGF) techniques, 
which are inherently quantum mechanical, but they suf- 
fer from the drawback of being suitable only to harmonic 
systems 4 ^—. 

Anharmonic systems on the other hand provide a tool 
to control heat transfer in nano systems with potential 
technological applications, e.g., a thermal diode£~— and a 
thermal transisto r 10 ' 11 . Some systems are purely anhar- 
monic like the spin boson model^, a paradigm in con- 
densed matter and quantum computation. Thus, it is im- 
portant to study the effect of anharmonic interactions in 
thermal transport from a fundamental point of view. To 
this end Wang developed quantum molecular dynamics 
which can probe into the moderate temperature regime 
but can not be extended to very low temperatures^. Se- 
gal et al. developed a master equation approach which 
is valid for all temperatures but their technique employs 
the Pauli master equatio n 12 ' 14 ' 15 which neglects possible 
coherent effects. Vclizhanin et al. combine the Green's 
function technique with the master equation approach 
but face the problem of non-conservation of energy 1 ^. 
Mingo 1 --, Wang et alJ^, and Santhosh. G et alJ£ have 
also developed techniques purely based on Green's func- 
tion which treat the anharmonicity perturbatively. 

Despite the various advances in molecular dynamics 
and NEGF techniques, the master equation approach 
seems the most suited tool to study thermal transport 
in anharmonic systems for arbitrary strength of anhar- 



monicity under weak coupling to the baths. In this ap- 
proach the heat current is calculated using the reduced 
density matrix along with an appropriate heat current 
operator. The reduced density matrix is typically cal- 
culated using a variety of master equations^— out of 
which the Redfield quantum master equation (RQME) 
is the most general equation with only the weak cou- 
pling approximation. The Pauli and the Lindblad master 
equations can be derived from it by introducing further 
approximations^. 

In this paper we will derive an explicit form of the heat 
current operator using standard perturbative techniques 
(for weak coupling) which along with the zeroth order 
contribution from the Redfield equation will allow us to 
calculate heat currents for anharmonic systems. One of 
the advantages of our formalism is that it can be used 
to study transients and it conserves energy in the steady 
state (without the need to symmetrize the heat current). 
In this work we will not address the problem of tran- 
sients; we will focus on the calculation of the steady state 
heat current for mono- and di-atomic molecules cither 
confined in a double well potential or connected by an 
anharmonic spring and having a non-linear on-site po- 
tential. 

For molecules confined in a double well potential we 
find that by varying the height of the barrier we ob- 
serve negative differential thermal conductance (NDTC) 
not only in the classical regime (high T) but also in the 
quantum regime (low T). However, in the problem with 
non-linear on-site potential and the anharmonic spring 
we find different behaviors of thermal conductance in 
all temperature ranges. Specifically at low-temperatures 
thermal conductance is drastically affected by the non- 
linear on-site potential whereas in the same temperature 
regime anharmonicity plays almost no role. 

The rest of the paper is organized as follows. In Sec.HIl 
we describe our basic model and the Redfield equation. 
In Sec. IIIII we derive the heat current operator using only 
the weak coupling approximation with standard quantum 
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mechanical perturbation theory. In Sec. IIVI wc discuss 
the bath and system models studied in this work. In 
Sec. [V] wc show our numerical results and comparisons 
with NEGF for the harmonic systems. Finally in Sec. IVII 
wc summarize our main conclusions. 



II. BASIC MODEL AND THE REDFIELD 
QUANTUM MASTER EQUATION 

Our basic model is similar to that used by many re- 
searchers for discussing thermal transport. It consists 
of a general system Hamiltonian connected to harmonic 
baths. The model Hamiltonian is thus of the Caldcira- 
Leggett type 2 ^, 
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where H s is some system Hamiltonian, a is a bath label 
allowing us to introduce multiple baths, P£ and Q% are 
the mass normalized normal variables of the a bath. The 
system-bath coupling strength parameter is e, and S a is 
the system operator coupled to the a th bath. In general 
it can be any function of the system operators. The above 
Hamiltonian can be rewritten as, 



where, 



1 k W * 

H" =S a ®B a . 



(2) 



Here B a = —tJ2k u kQk 1S ^ c collective bath operator 
that couples to the system. Throughout this paper we 
will set h — 1 and k B = 1. 

Assuming decoupled initial conditions for the total 
density matrix we can write the master equation for the 
reduced density matri x 24 ' 26 ' 27 as, 
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The relaxation coefficients are 

W$ = Wff + i (7 Q (0) + VF^J, (4) 

W* = J °dTC- lA -' T C a (r), (5) 

where 



Ay - Ei - Ej , 



C Q (r) = ^-r{m) (^coth f ^p) cos(wr) 

—i sui(u;t) I . (6) 



Ei is the i th energy of the system Hamiltonian and 
C(t) = (B a (r)i? a ) is the bath-bath correlator, where 
B a {r) is the free evolution operator according to 

exp[-ifl?r]. J Q H - ne 2 J2 k uf(2^)- 1 S(u J -^) 
is the spectral density of the bath, and 7 a (0) = 
7T _1 J dwJ Q (ct;)a; _1 is the damping kernel at time zero 
coming from the re-normalization part of the Hamilto- 
nian (H RN ). 

The above master equation is also known in the litera- 
ture as the Bloch-Rcdficld master equatio n 20 ' 24 ' 28 . While 
deriving it we have made only the weak coupling assump- 
tion. Other approximations such as the secular or ro- 
tating wave approximatio n 29 ' 30 or neglecting the Lamb 
shift a 31 ' 32 are commonly applied to Eq. ([3]). However 
these are uncontrolled approximations and we will not 
resort to them. 

We are primarily interested in studying the steady 
state heat current and hence to obtain the steady state 
reduced density matrix wc will set t — t = oo, and 
dp nm /dt = in Eq. ©. On the other hand the reduced 
density matrix in Eq. ([3]) can be formally written as a se- 
ries in the system-bath coupling parameter e, truncated 
at second order as, 



P = P {0) + 



eV 2 ). 



(7) 



Recently Fleming et alM- discussed that the density ma- 
trix obtained by setting dp nm /dt = in Eq. ([3]) gives in- 
accurate results for the second order diagonal elements. 
However, we will see in Scc. lIIII that the evaluation of heat 
current requires only p^ to which Fleming's argument 
does not apply. 

Now in order to obtain p^ we simply substitute the 
series Eq. ([7} in Eq. ^ and equate to zero the coefficients 
of all the powers of e. Thus by solving order by order we 
get the following equation for the steady state p(°> , 

E ( s n,sX - E szssfyA pL 0) = 0, (8) 

i,ot \ I / 



and for i ^ j we get, 



Pg>=0. 



(9) 
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Therefore using Eq. (J8j) and Eq. (|9|) along with the nor- 
malization condition Ti(p^ ) = 1 we can obtain the ze- 
roth order contribution to the reduced density matrix 



III. HEAT CURRENT 

Our system is connected with a semi-infinite left and 
a right heat bath^l whose Hamiltonians if£ and are 
defined in Eq. ([2) . In this section we will derive a formula 
for the heat current starting from the basic definition, 



dt 



(10) 



which is inspired by the change in energy of the (infi- 
nite) bath. The averaged operator is to be interpreted in 
the Heisenberg way, i dA/dt = [A, H tot ]. The time evolu- 
tion will be handled perturbativcly, much as one derives 
the RQME. Earlier works employing the master equation 
to calculate heat current have made additional approx- 
imations like symmetrization of the heat curren t 12 ' 16 or 
use of the Pauli master equation to calculate the reduced 
density matri x 12 ' 14 . Although this gives a simple form 
for the heat current operator, those approximations are 
not really needed, as will be shown here. 

Using the Heisenberg equation of motion in Eq. (|10[) 
we obtain, 



I L (i) 



(11) 



where 



A L (t) = (F L ®£ L )(i), 
F h = 5 L , 



(12) 



We recall that S L is the system operator connected to the 
bath operator B L of the left bath. The time evolution of 
the operator A L (t) is defined in terms of the evolution 
operator as, 



A^{t) = \J(t,t )^A L (t )V(t,t ), 



(13) 



Now we expand the evolution operator 
U(i, t a ) using the Kubo type identity^ 

order in e as^ 6 -, 



XJ(t,t 
U (*,t o 



e -iH (t-t ) 
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Here H" s (s) is the free evolution operator according to 
XJ (t,t o ). Using the above expression of the evolution 
operator in Eq. (fTB")) we get, 



A L (t) =A h (t)-i 



A^(t),H^(u) , (15) 



where A^(s) is again a free evolution. In order to obtain 
Eq. (fTS)) we have exploited the fact that the two heat 
baths are not directly coupled. We expanded only to 
first order since I L (t) in Eq. (fTTj) is already first order in 
e. 

From now on to simplify notation we will drop the bath 
label a. It is worth noting that even though Eq. (fT5|) has 
only the left bath label, the right bath comes in due to 
the free evolution of the operators ('.' H a = H s +J2 a H"). 
Now since in Eq. (|15p we require only the free evolution 
A(t) = F(t) <g) E(t) we express the operators F{t) and 
E{t) in terms of the free evolving Hubbard operator at 
time t as, 



(16) 



with X nm (t) = XJ (t,t o )l\m)(n\XJ (t,t o ), where \n),\m) 
are eigenvectors of the system Hamiltonian in the energy 
cigenbasis. Similarly, 

= E E Su9%n(u; t)X nm (t), (17) 



kl nm 



where 



g% m (u;t) = Tr [(X mn yui(u,t)X kl XJ (u,t)\ , (18) 

is a sort of freely evolving Green's function of the system. 

Now the operator A(t) can be expressed in terms of 
X(t) using Eq. (HU) and Eq. ([TTJ) in Eq. ([TSj) as, 

eA(t) = eY J X nm {t)®F nm E{t) 

( 



du 



X nj (t) ® F ru S kl g--E(t)B(u) 



\n,k,l 



- X lm {t) ® F jm S kl g- l 3 B(u)E{t) 



i>3 

n,k,l 



(19) 
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Using the factorized initial condition (p tot (t ) = p^(t ) > 
Ps(t a ) ® Pb(^o)) and tracing we obtain, 

e(A)=e^2(X nm (t))F nm (E(t)) 

n,m 

-i J2 (X nm (t)) J2 (Fn 3 S> m S<F jm ) , 



(20) 



where, 



— Sij 



due 



-i A k iu 



s < = (s>)\ 

X (u) = e 2 Tr B (E(t)B(t 



u)Pi 



(21) 
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and wc have used g^ m (u;t) — e l 5k, n ^l,m- Fi- 
nally, noting that (-B) = gives (E(t)) = 0, the heat 
current in Eq. (fTT|) (without the left bath label "L" ) can 
be expressed as, 



I = Tr(p(°)(t)l), 
1 = i (SS> - S<S) 



(22) 



where S< and S> are defined in Eq. ((2TJ) and p (a) (t) 
is the lowest order contribution to the reduced density 
matrb*2£. This is one of the main results of this paper. 
In order to evaluate the transients in heat current we 
require p^ and 1 at time t. I can be evaluated as long 
as we know the operators S > ' < at time t. p^°\t) can be 
calculated using the RQME (Eq. ©) by taking e very 
small while evaluating the bath-bath correlators C(r). 

Clearly from Eq. (|2"2"|) I = 2J and hence the 
heat current I is real. The correlator x( T ) 
:; f °° doju}J(u}) ^coth (j^j sin (am) + i cos(uju)j enter- 
ing in the current operator T can be expressed in terms 
of the bath-bath correlator C(r) (Eq. ([6])) used in the 
RQME since x(r) is the derivative of C(r). Therefore 
the operator S < can be computed as, 



s ij 



Sa C(0)-e 



-«Ajj({- to) 



C(t-t a ) 



iAijWij 



(23) 



Note that nothing particular to the harmonic baths has 
been invoked. Any other bath, e.g. spin baths, can be 
used as long as we can compute its bath correlators C(r) 
and x( T )- Thus only the relaxation rates W and the 
operators S < , S > are affected. 

Now using Eq. (f2"2"j) and Eq. (j2"3")l the heat current can 
be calculated in the steady state as well as in the transient 
where the S < operator has an explicit time dependence 
(some researchers refer to this as non-Markovian) . In this 
work we are interested in the steady state heat current 
and as mentioned in Sec. |TT] we will set t — t = oo. Since 
the bath correlator decays with time C(oo) will be zero 
for the steady state problem and thus only the transition 
rates W will contribute to the the current operator I (see 
Eq. {23])). Therefore using Eq. ®, Eq. ©, Eq. (p2|) . 
and Eq. (|2"3")l wc can calculate the heat current flowing 
through the system. 



IV. MODELS FOR SYSTEM AND BATHS 

A. Bath Model 

In order to describe the bath completely we need to 
specify a spectral density J(ui) which contains the in- 
formation about the frequency distributions of the bath. 
Several forms of the spectral density are used in the liter- 
ature based mainly on phcnomcnological modeling. Al- 
though the theory outlined in this work is not restricted 



to any particular form of the spectral density we will con- 
centrate on the Ohmic spectral density with a Lorentz- 
Drudc cut-off of the form, 



i + (w/cj D y 



(24) 



where r\ is the system bath coupling strength squared 
(e 2 ). One of the main advantages of using this spec- 
tral density is that we can calculate the bath correlators 
C(t) and the relaxation rates W analytically. Numerical 
decomposition of the spectral density is also used to an- 
alytically obtain the bath correlators, which reduces the 
computational cost a 38 i 39 . 

Decomposing the hyperbolic cotangent in Eq. into 
its Matsubara frequencies (yi = 2irlT) and noting that 
the resultant equation has poles at u = iiw D and 
w = ±i vi we can calculate the bath correlator using the 
theorem of residues as, 



C(t) = ^cot 



8 ^ 1 



vi e 



i=i 



-^e-^sginV). 



(f;K) 2 

(25) 



Once the bath correlator is obtained we can easily obtain 
the relaxation rates (W) using Eq. (j4j and Eq. ((5|) as, 



W>. = ^ 

13 2K + A?-) 



w D cot 



/3w D 



-A, 



13 2K+M.) 



B ^(i-(^ D ) 2 )(^ + A|.)' 

Bu \ uJu_ 
2 J Aij_ 



(26) 



cot 



7(0) = — ■ 



(27) 



(28) 



B. System 

Throughout our derivation of the formula for the heat 
current we have not specified the Hamiltonian of the sys- 
tem. In this work we will study the following systems; 
Model la: A monatomic molecule confined in a dou- 
ble well potential (Fig. (JTJl) ) , Model lb: A monatomic 
molecule with a linear + quartic non-linear on-site po- 
tential (Fig. (Ub)), Model 2a: A diatomic molecule con- 
fined in a double well potential where the atoms interact 
via a harmonic + quartic anharmonic spring (Fig. {It)), 
Model 2b: A diatomic molecule connected by a harmonic 
+ quartic anharmonic spring having a quartic non-linear 
on-site potential (Fig. {Th 1 )). 

In our models both on-site and coupling spring can 
be non-linear. However, we will restrict the use of non- 
linear to the on-site potential and use anharmonic for the 



couplings. The Hamiltonian of the monatomic molecule 
(Model la and lb) connected linearly via the position 
operator (S a — x; a = L, R) to two heat baths as shown 
in Fig. (QJi) and Fig. ((lb) is given by, 



if. 



+ \ a x A 



(29) 



where oj = \fk~o~ and we have set the mass of the atom 
to unity. When S = — 1 the potential has a double well 
structure (Duffing oscillator). When 5 = +1 (0 4 model) 
the system has a linear on-site potential with spring con- 
stant k plus a quartic non-linear term whose strength is 
governed by a parameter A . 

The Hamiltonian in case of the diatomic molecule 
(Model 2a and 2b) connected linearly via the position 
operator (5 L — x±, S R = X2) to two heat baths (Fig. (QJ) 
and Fig. (JTJi) ) is given by, 



i=l,2 
)2 



n 2 ( Xl - x 2 y 



+ X(xi - X2) 4 



(30) 



where w = ^/k~o, f2 = Vk and we have set the mass of the 
atom to unity. When 5 = — 1 similar to the monatomic 
case the diatomic molecule is confined in a double well 
potential and the atoms interact via a harmonic + an- 
harmonic spring. In case of S = +1 the model is gen- 
erally referred to as the FPU-/3 model and the atoms in 
the system are connected to each other via a harmonic 
spring k and a quartic anharmonic spring governed by a 
parameter A. Also each atom is subjected to a linear on- 
site potential whose spring constant is k and a quartic 
non-linear on-site potential whose strength is given by a 
parameter A . 

Model lb and Model 2b are of particular interest since 
they represent phonon-phonon interactions and to our 
knowledge these models have not been studied till date 
from the quantum (low temperature) to the classical 
regime (high temperature) for strong non-linearity or an- 
harmonicity. On the other hand Model la and Model 
2a are completely non-linear models and as we will see 
later exhibit interesting properties like negative differen- 
tial thermal conductance (NDTC) in quantum as well 
as classical regimes, which is a basic ingredient to build 
phononic devices like a thermal diode^. 
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FIG. 1. (Color Online) An illustration of the different systems 
connected to two heat baths at different temperatures T L (Red 
wave) and T R (Blue wave), (a) A monatomic molecule con- 
fined in a double well potential (Model la), (b) A monatomic 
molecule having a non-linear on-site potential (Model lb) . (c) 
A diatomic molecule confined in a double well potential and 
the atoms are connected by an anharmonic spring (Model 2a). 
(d) A similar diatomic molecule having a non-linear on-site 
potential (Model 2b). 



a system Hilbert space of around 40 levels is large enough 
to reach around 5 times the Debye temperature, i.e., suf- 
ficiently into the classical regime, whereas in case of the 
diatomic molecule a size of around 1600 (40x40) levels 
is sufficient to cover the same temperature range. 

In junction systems, since the cross-sectional area of 
the system interacting with the bath is not well defined, 
we can not define the thermal conductivity of the system. 
Hence in such cases we define thermal conductance as, 



lim 

t l ->t,t r -).t T l - T R 



(31) 



In order to numerically evaluate the thermal conductance 
we choose a small temperature difference between the two 
baths such that the limit in Eq. (f3"Tj) becomes valid. For 
all the systems considered in this work we find that a 
temperature difference of 10% is optimal and even if we 
decrease the temperature difference further the conduc- 
tance of the system does not change. 



V. RESULTS FOR THE HEAT CURRENT 



A. Heat current for the monatomic molecule 



1. Duffing oscillator model 



C. Numerical details 

In the numerical implementation, we can only use a 
finite number of base vectors. We will therefore choose 
a system Hilbert space large enough so that even at the 
highest temperatures the probability of finding the par- 
ticles in the highest energy levels is approximately zero. 
We do this by iteratively increasing the size of the system 
Hilbert space until at least five energy levels have a popu- 
lation less than 10 -15 . In case of the monatomic molecule 



We will first look at the Duffing oscillator 5 = — 1 
(Model la). Since it is a double well under no limiting 
case can this model be compared to exact solutions. Be- 
sides to the best of our knowledge it has not been looked 
at from the point of view of thermal transport. 

Fig. ([2]) shows the behavior of thermal conductance 
(Eq. ([31])) as a function of temperature (T = (T L +T R )/2) 
for varying strengths of the quartic term in the poten- 
tial A . First let us look at the two extreme cases when 
A /fc = 0.01 and \ /k a = 10 ((A 2 amu) _1 ). The barrier 
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Temperature (K) 

FIG. 2. (Color Online) Graph of conductance (a) vs temper- 
ature (T = (T L +Tr)/2) for various strengths of non-linearity 
(A ) in a monatomic molecule confined in a double well po- 
tential connected with Lorentz-Drude heat baths (Model la). 
Parameters used for the calculation are 8 — —1, k Q — 60.321, 
e = 6.0321 (meV/(A 2 amu)), uj d = 10 eV, and T R = 0.9T L . 
Ao/fco has dimensions of (A amu) -1 . 

height of the double well potential is inversely propor- 
tional to A and thus when \ a /k = 0.01 (A amu) -1 the 
particle remains confined to either one side of the bar- 
rier (indicated by the nearly degenerate eigenvalues in 

table (H|), whereas in case of X /k = 10.0 (A' amu 
the barrier is so low that the molecule simply experiences 
a quartic potential. Both these cases may be considered 
as the molecule experiencing only a quartic on-site po- 
tential. In these cases no NDTC behavior is observed in 
the quantum or classical regime. 

For intermediate values of the quartic term we observe 
negative differential thermal conductance (NDTC) be- 
havior in both the quantum and classical regimes. The 
main reason for this behavior is because for these values 
of the quartic strength the double well barrier is neither 
too strong nor too weak and hence the molecule can tun- 
nel through the barrier. At low temperatures and for 
certain intermediate values of quartic strength (X /k„ — 

0.05; 0.10 (A 2 amu) -1 in Fig ©) we see NDTC. In order 
to explain this possibly quantum behavior we analyze the 
lowest three eigenvalues and their populations given by 
as tabulated in table (JXJ> . Since the maxima of the 
double well potential barrier is at 0.0 eV clearly we can 
see from table dJ) that for X /k = 0.05; 0.10 (A 2 amu)" 1 
the lowest three eigenvalues are just below the maxima 
of the barrier indicating that the molecule can tunnel 
through the barrier and is not confined well within the 
double well as in the case of X a /k = 0.01 (A amu) -1 . 
The populations are also concentrated in the lowest two 
cigen states, supporting our claim. 

Now looking at the specific case of X /k = 

0.05 (A amu)" 1 we find that the lowest two energy lev- 



(a) (b) 




Left Bath Temperature (K) Left Bath Temperature (K) 

FIG. 3. (Color Online) Graph of current (I L ) vs tempera- 
ture of the left lead (T L ) using Landauer formula (Black) 
and our heat current formulation (Red) for the Lorentz- 
Drude model. The insets show current as a function of the 
strength of the dimensionless system-bath coupling squared, 
(a) shows the current comparison for a harmonic monatomic 
molecule (Model lb: S = +1, A = 0) and (b) shows the 
comparison for a harmonic diatomic molecule (Model 2b: 
S = +1, A = 0, A = 0). The parameters used for the 
monatomic molecule are; k = 60.321 meV/(A amu). The 
parameters used for the diatomic molecule are; k = 30.1605, 
k = 30.1605 (meV/(A amu)). The common bath parameters 
are; e = 6.0321 meV/(A 2 amu), uj d = 10 eV, and T R = 0.9T L . 
For both the insets the same system parameters and bath 
parameters are used except T L = 350 K and e 2 is varied. 



els are quite close, ~ 12.5meV (130 K). This is the ex- 
act temperature range at which the thermal conductance 
increases sharply indicating that the bath modes corre- 
sponding to that energy difference start conducting heat. 
In between 100 K to 300 K since the third energy level 
is quite far apart only the modes having energy corre- 
sponding to the energy difference of the first two energy 
levels transfer heat and thus the heat current nearly sat- 
urates showing NDTC. This claim is also supported by 
looking at the populations at 210 K which indicate that 
the third level has now started gaining some finite popu- 
lation. Then other energy modes will be allowed through 
the system causing the thermal conductance to again in- 
crease with temperature above 300 K. 

For other intermediate values of X /k = 

0.5; 1.0 (A' amu) the lowest lying energy states 
are above the potential barrier and hence in the quan- 
tum regime they do not experience the effect of the 
barrier. However, in the classical regime the molecule 
experiences the full potential, which includes the small 
barrier. Hence in the quantum regime no NDTC is 
observed, whereas a small NDTC is present in the 
classical regime. 
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Temperature (K) 

FIG. 4. (Color Online) Graph of conductance (a) vs temper- 
ature (T = (T L +Tr)/2) for various strengths of non-linearity 
in a monatomic molecule connected with Lorentz-Drude heat 
baths (Model lb). Parameters used for the calculation are 
8 = 1, k = 60.321, e = 6.0321 (meV/(A 2 amu)), uj d = 10 eV, 
and T R = 0.9T L . A /fc has dimensions of (A 2 amu) _1 . 

2. 4 model 

Now we will look at another model of the monatomic 
molecule known as the 4 model or the quartic on-site 
potential model (Model lb; 5 = +1). This model can be 
physically realized as a monatomic molecule interacting 
via non-linear interaction with a substrate. For this quar- 
tic on-site potential model in the limiting case of A = 0, 
the system becomes purely harmonic and we can employ 
NEGF techniques using the Landauer formula to evaluate 
the heat current as shown in Appendix. The Landauer 
formula is applicable only in the steady state for har- 
monic systems but its advantage is that it is applicable 
for all coupling strengths. Fig shows the heat current 
I L calculated via Landauer formula (Black curve) and our 
heat current formulation of Sec. HIT] (Red curve). The in- 
set shows the heat current as a function of the dimension- 
less system-bath coupling strength squared (e 2 /kg). We 
see that both curves exactly overlap in the weak system- 
bath coupling regime, i.e., up to e = 0.1fc o for the entire 
range of temperature showing excellent agreement be- 
tween the two approaches in this limit. For this specific 
harmonic case we have also compared our work to that of 
SegalM who obtains an analytic formula for heat current 
and we find excellent agreement between her approach 
and ours. 

Next we study the behavior of thermal conductance as 
a function of temperature (T = (T L + T R )/2) for varying 
strengths of non-linearity as shown in Fig. (Ql . It can be 
seen that even with the slightest amount of non-linearity 
the system behaves quite differently as compared to the 
harmonic case. The non-linearity does not change the 
behavior only in the high temperature (classical regime) 



TABLE I. Table of first three eigenvalues and correspond- 
ing populations for the mono- (A' = X a /k a ) and di-atomic 
(A' = A /(feo — k)) molecule confined in a double well poten- 
tial (Model la and Model 2a). 

A' Eigenvalues Populations in% 



Monatomic Diatomic Monatomic Diatomic 

(A 2 amu)" 1 (10~ 3 eV) (T = 210 K) (T = 105 K) 

-449.54 -1815.01 49.78 49.99 

0.01 -449.53 -1815.01 49.77 49.00 

-356.11 -1702.86 0.22 0.00 

-130.41 -424.28 67.33 50.25 

0.05 -117.94 -424.19 32.60 49.74 

-6.26 -287.59 O05 O00 

-74.89 -214.66 86.13 59.21 

0.10 -43.44 -211.45 13.84 40.77 

75.73 -85.70 O01 0.01 

10.60 18.01 99.53 99.96 

0.50 103.44 88.18 0.46 0.03 

242.11 135.98 O00 O00 

39.79 78.57 99.92 99.99 

1.00 163.29 186.79 0.07 0.00 

320.20 217.78 O00 O00 

133.16 266.28 99.99 99.99 

10.0 390.44 518.69 0.00 0.00 

688.89 528.42 0.00 0.00 



but also changes the behavior of low-temperature ther- 
mal conductance (quantum regime). To the best of our 
knowledge the effect of strong non-linearity in the quan- 
tum regime has not been studied and this simple system 
demonstrates that even in the low-temperature regime 
the non-linear forces can not be neglected. 



B. Conductance for the diatomic molecule 

Similar to the monatomic case we will first look at the 
case of the diatomic molecule trapped in a double well 
potential (6 — —1) where the atoms of the molecule inter- 
act only via a harmonic interaction (A = 0). We vary the 
height of the barrier by varying A and plot the conduc- 
tance as a function of temperature in Fig. <|5j) . Similar 
to the monatomic molecule case we observe NDTC in 
the quantum as well as classical regime for intermediate 
values of the strength of the non-linear potential. 

An analysis similar to the monatomic molecule case 
can be made with the eigenvalues and the pop- 
ulations shown in table |pj. For X /(k — k) = 

0.01,0.05 (A amu) -1 the barrier is very high and hence 
the molecule remains confined to either one side of 
the well indicated by the nearly degenerate eigenvalues 
and corresponding 50-50% probabilities (table (HJ). For 

A /(fco — k) = 0.1 (A amu) -1 we can observe NDTC in 
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FIG. 5. (Color Online) Graph of conductance (it) vs tem- 
perature (T — (T L + T R )/2) for a diatomic molecule con- 
fined in a double well potential using the Lorentz Drude bath 
model (Model 2a). The parameters used for the calcula- 
tion are: 5 = -1, k = 90.4815, k = 30.1605, e = 6.0321 
(meV/(A 2 amu)), A = meV/(A 4 amu 2 ), uj u = 10 eV, and 
Tr = 0.9T L . \ /(k — k) has dimensions of (A 2 amu) _1 . 

the quantum regime because only for this value the bar- 
rier is neither too high nor too low and hence the molecule 
can tunnel through the barrier since the eigenvalues are 
just below 0.0 eV (barrier maxima) as seen from table (jl|. 
The diatomic molecule brings another interesting aspect, 
i.e., the role of anharmonic interactions between the two 
connecting atoms. The anharmonic spring (A ^ 0) plays 
a small role in determining whether the system shows 
NDTC or not and it simply shifts the thermal conduc- 
tance in the high-temperature regime to a lower value as 
compared to the harmonic spring. This behavior is some- 
how expected since anharmonic interaction between the 
atoms leads to more scattering causing the thermal con- 
ductance to decrease as compared to a harmonic inter- 
action. Thus by looking at the monatomic and diatomic 
case it seems that the NDTC behavior in such double 
well potentials can be tuned by solely varying the barrier 
height (A ) and the number of atoms or anharmonicity 
in the system seem to play a small role. 

Now we will look at the FPU-/3 model with 5 = +1 
where we will first compare the heat current in a purely 
harmonic system (A = and A = 0) by our heat current 
formulation to the Landauer formula (Fig. ((Sb))- The 
inset shows the current as a function of the dimension- 
less system-bath coupling strength squared (e 2 / (k +k) 2 ). 
Again for weak system-bath coupling, i.e., up to e = 
0.1 (fc + k) there is an excellent agreement between our 
heat current formulation (Red) and the Landauer for- 
mula (Black) on the entire temperature range. 

Now we first switch on the anharmonicity, i.e., vary 
the parameter A and set the non-linear on-site poten- 
tial to zero i.e A = 0. Fig. © shows the effect of 
anharmonicity on the thermal conductance. Compar- 



FIG. 6. (Color Online) Graph of conductance (a) vs temper- 
ature (T = (T L +T H )/2) for an anharmonic diatomic molecule 
using the Lorentz Drude bath model (Model 2b). Parameters 
used for the calculation are: k — 30.1605, k = 30.1605, e = 
6.0321 (meV/(A 2 amu)), A = meV/(A 4 amu 2 ), uj d = 10 eV, 
and T R = 0.9T L . A/(fc + k) has dimensions of (A amu) -1 . 

ing Fig. ((4]) and Fig. ([6]) we see that the behavior of 
thermal conductance as a function of temperature for 
a non-linear on-site model and an anharmonic model 
is very different. For example in the low-temperature 
regime the behavior of thermal conductance for an an- 
harmonic diatomic molecule is same as that of a har- 
monic diatomic molecule. In Fig. ([7]) we plot the low- 
temperature thermal conductance for some combinations 
of non-linear on-site potential and anharmonicity for the 
diatomic molecule. Clearly only when we have non- 
linearity present the behavior of low-temperature ther- 
mal conductance differs from the harmonic case proving 
that the low-temperature behavior of thermal conduc- 
tance is strongly affected by non-linearity, whereas an- 
harmonicity plays almost no role at low temperatures. 

We would like to end this section with a few words on 
the specific heat of the systems considered in this work. 
Since specific heat of a solid is closely related to its ther- 
mal conductivity (k = 1/(3V) J2 k c kVjk, where c k is the 
specific heat, v k is the phonon group velocity, and l k is 
the mean free path associated with mode k and V is the 
volume of the solid) one expects a similar relation should 
hold even for the thermal conductance (a) . The quantum 
correction method 4 ^ - — is typically employed to relate the 
thermal conductance and the specific heat of the system. 
Although this approximation might not be valid for all 
temperatures 4 ^, in the low-temperature regime since the 
group velocity (v k ) can be approximated as a constant we 
get a oc C v . In case of the </> 4 model the low temperature 
specific heat^i shows similar behavior to the thermal con- 
ductance but at high temperatures the specific heat of ^ 4 
model shows a negative slope which is not observed in the 
thermal conductance indicating the role of the phonon 
group velocity. In case of mono- and di-atomic particle 
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FIG. 7. (Color Online) Graph of conductance (a) vs temper- 
ature (T = (T L + T R )/2) for an anharmonic + non-linear di- 
atomic molecule using the Lorentz Drude bath model (Model 
2b). Parameters used for the calculation are same as Fig. ([6]). 
<V(fco + k) and A /(fc + k) have dimensions of (A 2 amu) -1 . 

in a double well potential the physics behind the double- 
peak structure in the thermal conductance is similar to 
the Schottky anomaly of specific heat, which has been 
mainly studied for magnetic system o 47 ' 48 . Thus in case 
of anharmonic systems by observing the specific heat of 
materials it might be possible to predict features in the 
thermal conductance making it easy to choose materials 
for phononic devices. 



VI. CONCLUSION 

We have presented a fully quantum-mechanical "non- 
Markovian" theory based on standard perturbation the- 
ory to evaluate heat current in general anharmonic sys- 
tems. Our theory is valid for any strength of anhar- 
monicity and can be easily applied to any potential as 
long as the system- bath coupling is weak, i.e., up to 10% 
of the spring constant of the harmonic oscillator. Using 
this method we investigated thermal transport in mono- 
and di-atomic molecules confined in a double well po- 
tential and found that in this purely non-linear model 
we can tune the NDTC by simply varying the height of 
the barrier which is essential to make phononic devices 
to control the heat current. We also investigated the 
monatomic molecule having a non-linear on-site poten- 
tial and found that non-linearity affects the thermal con- 
ductance not only at high temperatures but also at low 
temperatures. In case of the diatomic molecule we found 
that in the low-temperature regime anharmonicity plays 
no vital role and the behavior of the thermal conductance 
is similar to the harmonic case. In order to quantify this 
statement we added a non-linear on-site potential to our 
diatomic molecule and found that the low-temperature 
thermal conductance deviated from the harmonic case 



proving that at low temperatures non-linearity can dras- 
tically affect the thermal conductance of the system. 

The technique presented here allows us to deal with 
any form of the system potential enabling us to study not 
only the phonon-phonon interactions from a fundamental 
point of view but also to explore systems of potential 
technological interest from the point of view of phononics. 
Although we have laid stress on the steady state thermal 
conductance in this work one important aspect in the 
field of phononics would be to study the effects of external 
fields on transient behaviors of purely non-linear systems 
which would enable us to control heat current and build 
better phononic devices. 
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APPENDIX: LANDAUER FORMULA AND 
MEAN-FIELD LIKE APPROXIMATION 

In case of harmonic systems the heat current can be 
obtained exactly for any arbitrary strength of the cou- 
pling using the Landauer formula^ given by, 

I L = J a 2^TM(/ L -/a), (32) 

where / LjR = (cxp[u;/T L |R ] — is the Bosc-Einstcin dis- 
tribution for phonons, and T[u)] is known as the transmis- 
sion coefficient. The transmission coefficient for a micro- 
scopic model is typically obtained by using the formula 
proposed by Meir et al.— given by, 

TM = Tr(GT L G<T R ), (33) 

where G r — (G a ) <! is the retarded Green's function 
and T L R describes interaction between the baths and 
the system. For a harmonic spring model described by 
the Caldeira-Leggett Hamiltonian (Eq. (JTJ) the retarded 
Green's function and r LR are given by, 



v ' ufl - K s - 2 ( 7l (0) + 7 r(0)) - Tr(u) ' 
T L , R = -2Im[E[ iR (w)], (34) 
S» = E[(u;) + E R (a;), (35) 

where K s is a spring constant matrix of the system having 
dimensions N x N, where N is the degrees of freedom of 
the system, 7 Q (0) is the damping kernel at time zero and 
£[ R (w) is the retarded self-energy of the left and right 
baths 4 ^. In case of a ID problem only one element of 
E[ R (w) matrix is non-zero and is given by, 

KM = -p / — [ —rdu;i - i J L '», (36) 



10 



where J L,R (w) is the spectral density of the bath. Using 
Eq. (|3"2"j) the thermal conductance defined in Eq. (j3"Tj) is 

given by, 



fd^ ,3/ 



(37) 



If the system Hamiltonian consists of a single harmonic 
oscillator (Eq. (|29|) ) and both the baths couple to the sys- 
tem with the same spectral density i.e J L (w) = J R (w) = 
J(oj) the transmission co-efficient is given by, 



TM 



4J 2 (w) 



(cj 2 - k ren - Re(E r (aj))Y + 4J 2 (w) 



(38) 



where fc rcn = k a + 47(0) is the re-normalized spring con- 
stant. Since we are interested in the weak coupling limit 
we make use of the following identity, 



lim 



. >6 (a; 2 — a 2 ) + e 2 
to obtain the transmission coefficient as 
7rJ(o; 



= ^-(6(x-a)+5(x + a)),(39) 



7>] 



(6(u - uj IBn ) + &(u + w te „)) , (40) 



where w ron = y/ fc ron . Therefore the thermal conductance 
in the weak coupling limit for a single harmonic particle 
in the system is given by, 



w ren e t 



rJ(w ro „). 



(41) 



One of the simplest approximations to treat non-linear 
on-site potential is a mean-field like approximation in 
which we transform the quartic problem into a quadratic 



one by replacing the quartic on-site potential in Eq. (|29[) 

by, 



(42) 



where 



<* 2 ) 



1 



a*. V coth (i?) +coth te )) ' (43) 



in the weak coupling limit. Here the particle in the 
quartic potential experiences a mean-field from its own 
quadratic part of the Hamiltonian making the system 
quadratic but the spring constant temperature depen- 
dent, 



2A 



coth 



( 2T ) 



(44) 



Although this is one of the crudest approximations it 
allows us to obtain the thermal conductance as, 



w mf e t 
2T 2 (e-r 1 -1 



J (W»f) 



(45) 



where w mf = \J k m! . In case of the harmonic system and 
the non-linear on-site model under the mean-field like 
approximation we can obtain the low-temperature de- 
pendence of the thermal conductance from Eq. (|4ip and 
Eq. (j45]) as, 



e t 



er wc cx 



Cwc mf °c 



rp2 ' 

y/fcren + (2* /\/' t rcn) 



J^2 



(46) 
(47) 



which clearly shows that even for the mean-field like ap- 
proximation the low-temperature thermal conductance 
for a quartic on-site model is quite different from the 
harmonic system. In comparison with our heat current 
formulation this mean-field like approximation gives the 
correct qualitative features for the thermal conductance 
although the exact behavior is different. 
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